Introducing the R Spatial Ecosystem

The basics

A few main packages are dedicated to spatial data handling like import/export, display, geoprocessing, spatial ststistics.

“Historical” Packages

rgdal: interface between R and GDAL (Geospatial Data Abstraction Library) and PROJ4 libraries.
sp: classes and methods for spatial data in R.
rgeos: interface between GEOS (Geometry Engine - Open Source) library and R: area, perimeter, distances, dissolve, buffer, overlap, union, contains…

These packages are still widely used.

Simple Features for R

sp, rgeos and rgdal functionnalities in one package.
Easier data handling, simpler objects.
Main author and maintainer: Edzer Pebesma (also sp author)
Compatibility with the pipe synthax and tidyverse operators.
First release: October 20, 2016
Website: Simple Features for R

sf Objects Data Structure

format sf

Using `sf`

Data Import

library(sf)
mtq <- st_read("data/mtq/martinique.shp")

Reading layer `martinique' from data source `/home/tim/Documents/prz/flo/data/mtq/martinique.shp' using driver `ESRI Shapefile'
Simple feature collection with 34 features and 23 fields
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: 690574.4 ymin: 1592426 xmax: 736126.5 ymax: 1645660
epsg (SRID):    32620
proj4string:    +proj=utm +zone=20 +datum=WGS84 +units=m +no_defs

Data Display

Default

plot(mtq)

Only geometry

plot(st_geometry(mtq))

Centroids Extraction

mtq_c <- st_centroid(mtq)
plot(st_geometry(mtq))
plot(st_geometry(mtq_c), add = TRUE, cex = 1.2, col = "red", pch = 20)

Distance Matrix

mat <- st_distance(x = mtq_c, y = mtq_c)
mat[1:5, 1:5]

Units: m
          [,1]     [,2]      [,3]      [,4]      [,5]
[1,]     0.000 35297.56  3091.501 12131.617 17136.310
[2,] 35297.557     0.00 38332.602 25518.913 18605.249
[3,]  3091.501 38332.60     0.000 15094.702 20226.198
[4,] 12131.617 25518.91 15094.702     0.000  7177.011
[5,] 17136.310 18605.25 20226.198  7177.011     0.000

Polygons Aggregation

mtq_u <- st_union(mtq)
plot(st_geometry(mtq), col = "lightblue")
plot(st_geometry(mtq_u), add = T, lwd = 2, border = "red")

Buffer Zone

mtq_b <- st_buffer(x = mtq_u, dist = 5000)
plot(st_geometry(mtq), col = "lightblue")
plot(st_geometry(mtq_u), add = T, lwd = 2)
plot(st_geometry(mtq_b), add = T, lwd = 2, border = "red")

Polygon Intersection

m <- rbind(c(700015, 1624212), c(700015, 1641586), c(719127, 1641586), c(719127, 
    1624212), c(700015, 1624212))
p <- st_sf(st_sfc(st_polygon(list(m))), crs = st_crs(mtq))
plot(st_geometry(mtq))
plot(p, border = "red", lwd = 2, add = T)

mtq_z <- st_intersection(x = mtq, y = p)
plot(st_geometry(mtq))
plot(st_geometry(mtq_z), col = "red", border = "green", add = T)

Voronoi Polygons

google: “st_voronoi R sf” (https://github.com/r-spatial/sf/issues/474 & https://stackoverflow.com/questions/45719790/create-voronoi-polygon-with-simple-feature-in-r)

mtq_v <- st_voronoi(x = st_union(mtq_c))
mtq_v <- st_intersection(st_cast(mtq_v), st_union(mtq))
mtq_v <- st_join(x = st_sf(mtq_v), y = mtq_c, join = st_intersects)
mtq_v <- st_cast(mtq_v, "MULTIPOLYGON")
plot(st_geometry(mtq_v), col = "lightblue")

Mapping

Static Maps

Several solutions are available:

ggplot2 users can have a look to ggplot2 mapping features (geom_sf) that can mix nicely with ggspatial.
For more advanced mapping features in a ggplot2-like syntax have a look to tmap
cartography is based on base graphics and allow most of basic and advanced cartographic representations.
Here we will focus on cartography. Full disclosure: the speaker is the maintener of cartography.

Using `cartography`

Data Preparation

library(sf)
# Import geo layer
sm <- st_read(dsn = "data/rhone/seine_maritime.geojson", stringsAsFactors = F)

Reading layer `OGRGeoJSON' from data source `/home/tim/Documents/prz/flo/data/rhone/seine_maritime.geojson' using driver `GeoJSON'
Simple feature collection with 718 features and 1 field
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: 0.06560943 ymin: 49.25139 xmax: 1.790225 ymax: 50.07121
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs

dep <- st_read(dsn = "data/rhone/dep.geojson", stringsAsFactors = F)

Reading layer `OGRGeoJSON' from data source `/home/tim/Documents/prz/flo/data/rhone/dep.geojson' using driver `GeoJSON'
Simple feature collection with 96 features and 2 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: -5.138996 ymin: 41.36216 xmax: 9.559598 ymax: 51.089
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs

# change projection (lambert93)
sm <- st_transform(sm, 2154)
dep <- st_transform(dep, 2154)

# Import dataset
csp <- read.csv("data/rhone/data_seine_maritime.csv")
# merge geolayer and dataset
sm <- merge(sm, csp, by = "INSEE_COM", all.x = TRUE)

Proportionnal symbols

library(cartography)
# Map of active population
plot(st_geometry(sm))
propSymbolsLayer(sm, var = "act")
title("Active Population")

# Custom map of active population
par(mar = c(0.2, 0.2, 1.4, 0.2), bg = "azure")
bb <- st_bbox(sm)
plot(st_geometry(dep), col = "ivory", border = "ivory3", xlim = bb[c(1, 3)], 
    ylim = bb[c(2, 4)])
plot(st_geometry(sm), col = "cornsilk2", border = NA, lwd = 0.5, add = T)
propSymbolsLayer(sm, var = "act", col = "darkblue", inches = 0.5, border = "white", 
    lwd = 0.7, symbols = "square", legend.style = "e", legend.pos = "topleft", 
    legend.title.txt = "Labor Force\n(2014)", legend.values.rnd = 0)
barscale(size = 10)
north(pos = "topright", col = "darkblue")
layoutLayer(title = "Workforce in Seine-Maritime", sources = "Insee, 2018", 
    author = "Kim & Tim, 2018", col = "darkblue", coltitle = "white", frame = TRUE, 
    scale = NULL, north = FALSE)

Typology Map

par(mar = c(0.2, 0.2, 1.4, 0.2), bg = "azure")
plot(st_geometry(dep), col = "ivory", border = "ivory3", xlim = bb[c(1, 3)], 
    ylim = bb[c(2, 4)])
typoLayer(sm, var = "cat", border = "ivory", lwd = 0.5, legend.values.order = c("agr", 
    "art", "cad", "int", "emp", "ouv"), col = c("#e3b4a2", "#a2d5d6", "#debbd4", 
    "#b5dab6", "#afc2e3", "#e9e2c1"), add = TRUE, legend.pos = "n")

legendTypo(title.txt = "Dominant Category", col = c("#e3b4a2", "#a2d5d6", "#debbd4", 
    "#b5dab6", "#afc2e3", "#e9e2c1"), categ = c("Agriculteurs", "Artisans", 
    "Cadres", "Prof. Inter.", "Employés", "Ouvriers"), nodata = F)
barscale(size = 10)
north(pos = "topright", col = "darkblue")
layoutLayer(title = "Workforce Distribution in Seine-Maritime", sources = "Insee, 2018", 
    author = "Kim & Tim, 2018", col = "darkblue", coltitle = "white", frame = TRUE, 
    scale = NULL, north = FALSE)

Choropleth Map

# Cartes choroplèthes
sm$pcad <- 100 * sm$cad/sm$act

bks <- getBreaks(v = sm$pcad, method = "quantile", nclass = 6)
cols <- carto.pal("green.pal", 3, "wine.pal", 3)
par(mar = c(0.2, 0.2, 1.4, 0.2), bg = "azure")
plot(st_geometry(dep), col = "ivory", border = "ivory3", xlim = bb[c(1, 3)], 
    ylim = bb[c(2, 4)])
choroLayer(sm, var = "pcad", breaks = bks, col = cols, border = "grey80", lwd = 0.4, 
    legend.pos = "topleft", legend.title.txt = "Part des cadres\n(en %)", add = T)
layoutLayer(title = "Les cadres", sources = "Insee, 2018", author = "Kim & Tim, 2018", 
    theme = "green.pal", col = "darkred", coltitle = "white", postitle = "center", 
    frame = TRUE, scale = 10)
north(pos = "topright", south = TRUE)

library(cartogram)
sm_c1 <- cartogram_cont(sm, "act", 3)
sm$pcad <- 100 * sm$cad/sm$act

bks <- getBreaks(v = sm$pcad, method = "quantile", nclass = 6)
cols <- carto.pal("green.pal", 3, "wine.pal", 3)
par(mar = c(0.2, 0.2, 1.4, 0.2), bg = "azure")
plot(st_geometry(dep), col = "ivory", border = "ivory3", xlim = bb[c(1, 3)], 
    ylim = bb[c(2, 4)])
choroLayer(sm_c1, var = "pcad", breaks = bks, col = cols, border = "grey80", 
    lwd = 0.4, legend.pos = "topleft", legend.title.txt = "Part des cadres\n(en %)", 
    add = T)
layoutLayer(title = "Les cadres", sources = "Insee, 2018", author = "Kim & Tim, 2018", 
    theme = "green.pal", col = "darkred", coltitle = "white", postitle = "center", 
    frame = TRUE, scale = 10)
north(pos = "topright", south = TRUE)

Gridded Map

These maps may allow to hide or, at least, diminish the MAUP.

sqrt(quantile(st_area(sm), 0.9))

3904.689 m

# Grid map
grid <- getGridLayer(x = sm, cellsize = 3904 * 3904, type = "regular", var = c("cad", 
    "act"))
grid$pcad <- 100 * grid$cad/grid$act
par(mar = c(0.2, 0.2, 1.4, 0.2), bg = "azure")
plot(st_geometry(dep), col = "ivory", border = "ivory3", xlim = bb[c(1, 3)], 
    ylim = bb[c(2, 4)])
choroLayer(grid, var = "pcad", breaks = bks, col = cols, border = "grey80", 
    lwd = 0.4, legend.pos = "topleft", legend.title.txt = "Part des cadres\n(en %)", 
    add = T)
layoutLayer(title = "Les cadres", sources = "Insee, 2018", author = "Kim & Tim, 2018", 
    theme = "green.pal", col = "darkred", coltitle = "white", postitle = "center", 
    frame = TRUE, scale = 10)
north(pos = "topright", south = TRUE)

Smoothed Map

# Smooth Map
grid$cad100 <- grid$cad * 100
par(mar = c(0.2, 0.2, 1.4, 0.2), bg = "azure")
plot(st_geometry(dep), col = "ivory", border = "ivory3", xlim = bb[c(1, 3)], 
    ylim = bb[c(2, 4)])
smoothLayer(x = grid, var = "cad100", var2 = "act", typefct = "exponential", 
    span = 3900, beta = 2, breaks = bks, col = cols, legend.pos = "topleft", 
    mask = st_buffer(sm, 0), legend.title.txt = "Part des cadres\n(en %)*", 
    border = "grey90", lwd = 0.2, add = T)
layoutLayer(title = "Les cadres", sources = "Insee, 2018", author = "Kim & Tim, 2018", 
    theme = "green.pal", col = "darkred", coltitle = "white", postitle = "center", 
    frame = TRUE, scale = 10)
north(pos = "topright", south = TRUE)
text(x = 518000, y = 6921123, font = 3, cex = 0.8, labels = "*Lissage par potentiels\n fonction exponnentielle\n span = 2.5km, beta = 2")

Dynamic/interactive maps

The 2 main packages dedicated to interactive mapping are leaflet and mapview.

leaflet

mapview

Spatial modelling and MISC

CRAN Task View: Analysis of Spatial Data

Exploration of sirene & osm db for restaurants

See file data_prep.R for data extraction.

# spatial data management
library(sf)
# cartography
library(cartography)
# smooth density computation
library(spatstat)
library(raster)
library(maptools)
# interactive mapping
library(mapview)

Display points

# load data
adm <- readRDS("data/iris_p13.rds")
feat_sir <- readRDS("data/sir_p13.rds")
feat_osm <- readRDS("data/osm_p13.rds")

nrow(feat_sir)

[1] 1066

nrow(feat_osm)

[1] 585

par(mar = c(0, 0, 1.2, 0), mfrow = c(1, 2))

plot(st_geometry(adm), col = "ivory1", border = "ivory3", lwd = 0.5, bg = "#FBEDDA")
plot(st_geometry(feat_sir), col = "red", pch = 20, add = T, cex = 0.5)
layoutLayer(title = "SIR", sources = "", author = "", scale = NULL, frame = FALSE)
legend(x = "topright", legend = "= 1 restaurant", pch = 20, pt.cex = 0.5, cex = 0.7, 
    col = "red")

plot(st_geometry(adm), col = "ivory1", border = "ivory3", lwd = 0.5, bg = "#FBEDDA")
plot(st_geometry(feat_osm), col = "red", pch = 20, add = T, cex = 0.5)
layoutLayer(title = "OSM", sources = "", author = "", scale = 0.5, north = T, 
    frame = FALSE)

Count points in admin units

inter_osm <- st_intersects(adm, feat_osm)
inter_sir <- st_intersects(adm, feat_sir)

adm <- st_sf(adm[, 1:2, drop = T], n_osm = sapply(X = inter_osm, FUN = length), 
    n_sir = sapply(X = inter_sir, FUN = length), geometry = st_geometry(adm))

par(mar = c(0, 0, 1.2, 0), mfrow = c(1, 2))
plot(st_geometry(adm), col = "ivory1", border = "ivory3", lwd = 0.5, bg = "#FBEDDA")
propSymbolsLayer(adm, var = "n_sir", inches = 0.1)
layoutLayer(title = "SIR", sources = "", author = "", scale = NULL, frame = FALSE)

plot(st_geometry(adm), col = "ivory1", border = "ivory3", lwd = 0.5, bg = "#FBEDDA")
propSymbolsLayer(adm, var = "n_osm", inches = 0.1)
layoutLayer(title = "OSM", sources = "", author = "", scale = 0.5, north = T, 
    frame = FALSE)

# plot(adm$n_sir, adm$n_osm, asp = T, xlim = c(0,55), pch = 21, bg = 'red')
# text(adm$n_sir, adm$n_osm, labels = adm$CODE_IRIS, cex = 0.5, pos = 2) x
# <-lm(formula = adm$n_osm~adm$n_sir) abline(a = x$coefficients )

Interactive visualisation

Explore datasets interactively

The default is not bad:

mapview(feat_sir) + mapview(feat_osm, col.regions = "red")

mapview is hightly customizable:

mapview(feat_sir, map.types = "OpenStreetMap", col.regions = "#940000", label = paste(feat_sir$L2_NORMALISEE, 
    feat_sir$NOMEN_LONG, sep = " - "), color = "white", legend = TRUE, layer.name = "SIRENE", 
    homebutton = FALSE, lwd = 0.5) + mapview(feat_osm, col.regions = "#000094", 
    color = "white", legend = TRUE, label = feat_osm$name, lwd = 0.5, layer.name = "OSM", 
    homebutton = FALSE)

count points in regular grid

fufun2 <- function(feat, adm, cs = 500) {
    grid <- st_make_grid(x = adm, cellsize = cs, what = "polygons")
    x <- st_intersects(grid, feat)
    gg <- st_sf(n = sapply(X = x, FUN = length), grid)
    plot(adm$geometry, col = "ivory1", border = "ivory3", lwd = 0.5, bg = "#FBEDDA")
    choroLayer(gg[gg$n > 0, ], var = "n", border = NA, breaks = c(1, 4, 6, 9, 
        12, 14, 17, 19, 22), col = carto.pal("wine.pal", 8), add = TRUE)
}

par(mar = c(0, 0, 1.2, 0), mfrow = c(2, 3))
fufun2(feat_sir, adm, 50)
fufun2(feat_sir, adm, 100)
fufun2(feat_sir, adm, 150)
fufun2(feat_osm, adm, 50)
fufun2(feat_osm, adm, 100)
fufun2(feat_osm, adm, 150)

smooth density analysis

Define a function that compute smoothed restaurant density and map it.

fufun <- function(feat, adm, title, sigma = 100, res = 50) {
    bb <- as(feat, "Spatial")
    bbowin <- as.owin(as(adm, "Spatial"))
    pts <- coordinates(bb)
    p <- ppp(pts[, 1], pts[, 2], window = bbowin)
    ds <- density.ppp(p, sigma = sigma, eps = res)
    rasdens <- raster(ds) * 1000 * 1000
    proj4string(rasdens) <- "+init=epsg:2154"
    rasdens <- rasdens + 1
    bks <- getBreaks(values(rasdens), nclass = 12, method = "equal")
    cols <- (mapview::mapviewGetOption("raster.palette"))(length(bks) - 1)
    plot(adm$geometry, col = NA, border = NA, main = "", bg = "#FBEDDA")
    plot(rasdens, breaks = bks, col = cols, add = T, legend = F)
    legendChoro(pos = "topright", cex = 0.7, title.cex = 0.7, title.txt = "resto per km2)", 
        breaks = bks - 1, nodata = FALSE, values.rnd = 0, col = cols)
    layoutLayer(title = paste0(title), scale = 0.5, tabtitle = TRUE, frame = FALSE, 
        author = "Map data © OpenStreetMap contributors, under CC BY SA. SIRENE 2018", 
        sources = "cartography 2.0.2")
    return(invisible(rasdens))
}

par(mar = c(0, 0, 1.2, 0), mfrow = c(2, 3))
fufun(feat_sir, adm, "SIR 50", 50, 20)
fufun(feat_sir, adm, "SIR 100", 100, 20)
fufun(feat_sir, adm, "SIR 250", 250, 20)
fufun(feat_osm, adm, "OSM 50", 50, 20)
fufun(feat_osm, adm, "OSM 100", 100, 20)
fufun(feat_osm, adm, "OSM 250", 250, 20)

reproducibility

sessionInfo()

R version 3.5.1 (2018-07-02)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Debian GNU/Linux 9 (stretch)

Matrix products: default
BLAS: /usr/lib/libblas/libblas.so.3.7.0
LAPACK: /usr/lib/lapack/liblapack.so.3.7.0

locale:
 [1] LC_CTYPE=fr_FR.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=fr_FR.UTF-8        LC_COLLATE=fr_FR.UTF-8    
 [5] LC_MONETARY=fr_FR.UTF-8    LC_MESSAGES=fr_FR.UTF-8   
 [7] LC_PAPER=fr_FR.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=fr_FR.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] mapview_2.4.0         maptools_0.9-2        raster_2.6-7         
 [4] sp_1.3-1              spatstat_1.56-0       rpart_4.1-13         
 [7] nlme_3.1-137          spatstat.data_1.3-1   SpatialPosition_1.2.0
[10] cartogram_0.1.0       cartography_2.1.1     sf_0.6-3             
[13] rmdformats_0.3.3      knitr_1.20           

loaded via a namespace (and not attached):
 [1] Rcpp_0.12.18         lattice_0.20-35      deldir_0.1-15       
 [4] png_0.1-7            class_7.3-14         rprojroot_1.3-2     
 [7] digest_0.6.15        foreach_1.4.4        mime_0.5            
[10] plyr_1.8.4           R6_2.2.2             backports_1.1.2     
[13] stats4_3.5.1         evaluate_0.11        e1071_1.6-8         
[16] tensor_1.5           highr_0.7            rstudioapi_0.7      
[19] miniUI_0.1.1.1       R.utils_2.6.0        R.oo_1.22.0         
[22] Matrix_1.2-14        goftest_1.1-1        rmarkdown_1.10      
[25] webshot_0.5.0        rgdal_1.3-3          stringr_1.3.1       
[28] questionr_0.6.3      foreign_0.8-70       htmlwidgets_1.2     
[31] munsell_0.5.0        polyclip_1.9-0       shiny_1.1.0         
[34] compiler_3.5.1       httpuv_1.4.4.2       xfun_0.3            
[37] base64enc_0.1-3      mgcv_1.8-24          rgeos_0.3-28        
[40] htmltools_0.3.6      bookdown_0.7         codetools_0.2-15    
[43] viridisLite_0.3.0    later_0.7.3          R.methodsS3_1.7.1   
[46] grid_3.5.1           jsonlite_1.5         spData_0.2.9.0      
[49] satellite_1.0.1      xtable_1.8-2         DBI_1.0.0           
[52] magrittr_1.5         formatR_1.5          scales_0.5.0        
[55] units_0.6-0          stringi_1.2.4        promises_1.0.1      
[58] leaflet_2.0.1        spatstat.utils_1.8-2 RColorBrewer_1.1-2  
[61] iterators_1.0.10     tools_3.5.1          gdalUtils_2.0.1.14  
[64] crosstalk_1.0.0      abind_1.4-5          yaml_2.2.0          
[67] colorspace_1.3-2     classInt_0.2-3

Geovisualization and Spatial Analysis